Integrated study of quaternary aquifer for hydrostratigraphy and groundwater quality assessment in central Thal Doab, Punjab, Pakistan

The groundwater resources in different areas of Pakistan are heading towards depletion along with the deterioration of quality due to over-abstraction and urbanization. The main focus of this study is to map the current hydrostratigraphical and hydraulic conditions of the late Quaternary aquifers in the central part of Thal Doab of Punjab Plains. To achieve the target, a comprehensive approach was employed combining geophysical investigations using electrical resistivity surveys (ERS) and physiochemical analysis of groundwater specimens collected from the study area. Careful calibration of resistivity models was performed by comparing them with lithologs to ensure their accuracy. The current groundwater conditions were assessed through thirty vertical electrical soundings (VES) using the Schlumberger electrode configuration up to 300m of AB/2. The interpreted results revealed the presence of four to six geo-electric sublayers comprising the intermixing layers of clay, silt, sand, gravel, and kankar inclusions. These layers exhibited very low (<20 Ω-m) to very high (>230 Ω-m) resistivity zones at various depth intervals. The developed 2D/3D models of aquifer systems identify the promising areas of good/fresh quality groundwater in the regions characterized by medium to very high resistivity mainly within the sand with gravel layers. However, lower resistivity values indicate the presence of marginally suitable/fair and saline/brackish groundwater showing the existence of fine sediments such as clays/silts. Additionally, twenty groundwater samples were collected to assess various parameters including pH, TDS, arsenic, fluoride, iron, nitrate, and nitrite. The spatial distribution of these parameters was visualized using 2D maps. The suitability of the groundwater for drinking consumption was evaluated in accordance with WHO guidelines.


Introduction
The Late Quaternary aquifer system has been prone to adverse changes in Punjab region for the last two decades due to groundwater depletion, water table declination, excessive groundwater mining, sporadic rainfall patterns, and especially due to the shortage of river flow across the Punjab [1][2][3][4][5].The aforementioned concerns are also being observed in the study area and its vicinities as well.The study area is the part of Punjab Plateform which holds a thick subsurface alluvial cushion.Furthermore, the study area is the western part of district Muzaffargarh which lies in the central part of Thal Doab (region between Chenab and Indus Rivers), Punjab province, Pakistan.The present research work is to investigate the latest groundwater scenario according to the subsurface hydrostratigraphy, hydrogeological conditions, and groundwater quality assessment of the study area.
Subsurface aquifers hold the majority of freshwater resources, highlighting the necessity to comprehend their hydrological characteristics to effectively manage groundwater resources and ensure a sustainable freshwater supply for the community [6][7][8][9].Accurate evaluation of lithological and geomorphological interconnections in sedimentary deposits is crucial for monitoring subsurface and groundwater resources effectively.However, defining lithological composition in alluvial environments with diverse lithology and complex geomorphic features poses significant challenges.Direct observational techniques like drilling wellbores are often impractical due to time, cost, and environmental constraints.Surface geophysical methods offer a more effective approach for investigating large spatial and temporal regions.These methods have found extensive application in numerous fields to assess the hydrogeological aspects of confined or unconfined aquifers [8,[10][11][12][13].The advancements in electronic technology and the development of digital simulation solutions for hydrogeological challenges have significantly enhanced the effectiveness of geophysical research methods for aquifer exploration and groundwater modeling in the recent past [14][15][16][17][18].
Notably, the electrical resistivity method is believed as one of the best geophysical methods in groundwater exploration.It is commonly utilized for groundwater monitoring due to its affordability, ease of implementation, and high efficiency particularly in regions of heterogeneous lithological units [17,19,20].This method is frequently employed to address a range of groundwater issues including locating the saline-fresh groundwater zones and their boundaries [21][22][23] and demarcation of groundwater contaminant zones [24][25][26].Due to the advancements in computational modeling techniques, this method has gained more popularity in recent years [27][28][29].Factors such as groundwater salinity, aquifer pore fluids, rock type, clay content, porosity, and mineralogical composition influence the electrical resistivity of alluvium [30][31][32].The lithology, geological structures, and porosity of subsurface material significantly impact the distribution of pore-water below the water table [3,33,34].The porosity of sand and silts is affected by sorting levels, grain morphology, and compaction [35].Clay, due to its high capacity for retaining charged particles and moisture, exhibits low resistivity values.
Previous research works have utilized various geophysical techniques to characterize different aquifer types.This study aims to characterize subsurface lithological units down to 300m of AB/2, providing groundwater resource zones and delineating the saline-freshwater boundaries.The research includes a comprehensive assessment integrating hydrogeochemical analysis, field inspections, and laboratory testing of groundwater samples to infer the current groundwater scenario.The findings of this study are crucial for the sustainable management of groundwater resources in the research region.A total of 30 sparsely distributed electrical resistivity survey profiles were acquired in the research region.Moreover, 20 groundwater samples were also collected to analyse the amount of essential physiochemical parameters.

Geology, hydrostratigraphy, and geomorphology of the study area
Punjab province is mainly located in Doabs, characterized by interfluves systems bounded by rivers.These Doabs are interconnected with the extensive network of the Indus River and its tributaries spanning approximately 1200 km in length [2,[36][37][38].The specific study area is situated in the central region of the Thal Doab and its surface is primarily comprised of floodplain deposits of the lower terrace and loess deposits of the upper terrace (Fig 2) [2,39].Thal Doab is flanked by the Chenab and Indus Rivers (Fig 2 ) which serve as the primary sources of groundwater recharge while the rainfall plays a secondary role.The deposition of sediments in the Thal Doab area resulted from the dynamic processes of current and ancient water channels originating from the adjacent Indus and Chenab Rivers.
The shallow geological units in Thal Doab are Quaternary alluvium constituting the subsurface aquifer (Fig 2).These alluvial deposits predominantly consist of loosely consolidated fine to coarse sand, occasionally intercalated with discrete lenses of intermixed fine materials such as silts and clays.These sediments were deposited by the Indus River and its tributaries including Jhelum and Chenab Rivers within a subsiding depression.Thal Doab also exhibits extensive surface coverage of aeolian sand which overlays the alluvium.The thickness of alluvium varies as a thinner cushion in the central part of Thal Doab and a thicker column in peripheral sections adjacent to the mountainous areas.This alluvium overlays sedimentary strata which pinch out towards the east and increase towards the west and southwest.Furthermore, the sedimentary strata overlie the Precambrian Basement rocks of igneous/metamorphic nature.The northern boundary of Thal Doab is marked by the Salt Range, characterized by highly deformed and fossiliferous rock formations ranging from Pre-Cambrian to Pleistocene ages [39].Sedimentary rocks extending from Kirther and Suleman Ranges, border the western edge of Thal Doab [37,39-41,].The borehole data suggests that a >300m thick cover of alluvium comprised of intermixed layers of fine-grained silt/clay to medium to coarse sand with Late Quaternary gravel is present throughout this area.The sediment cover is heterogeneous with significant lateral and horizontal variation in lithology [42].Some of the borehole lithological logs representing the subsurface geology of the study area are shown in (Fig 3).
The aquifer in Thal Doab is composed of intermixed layers of coarser sediments like sand and gravel representing the primary form of available groundwater storage.The hydrogeological characteristics significantly influence the supply and replenishment of the aquifer.Finer deposits such as clayey and silty contents found in the subsurface can be presumed as a result of reworked loess deposits [2,39,41].Both the coarser and finer materials are porous sediments, however, the permeability of coarser ones such as sand and gravel exceed that of fine ones.Thus, the fine content acts as an aquitard, impeding water flow and contributing to natural groundwater salinity [42,43].Various parts of Thal Doab face hindrances to groundwater flow due to the intermixing of such fine sediments.The Indus River with its coarse grained sediments and higher elevation plays a pivotal role in regulating hydrology and drainage in Upper Thal Doab, facilitating the aquifer replenishment.The drainage system fed by numerous off-takes of the Indus River traverses Thal Doab which is a recharging source of alluvium aquifer through a natural gradient.Deterioration of groundwater quality can also be observed; however, its quality tends to improve in closer proximity to the river beds [2,4].A shallow water table prevails across most of the area typically ranging from 0.5 to 9 m.The digital elevation map and variations in groundwater levels across the study area showed a correlation with a decline in elevation from northwest to southeast (Fig 4).The lithological units of Indus River deposits exhibit high transmissivity leading to accelerating the groundwater recharge.

Geophysical investigation techniques
Over time, the electrical resistivity survey (ERS) has been widely utilized for decades as an effective tool for detecting and characterizing subsurface hydrogeological features [44][45][46][47].Apart from providing information about the lithological composition of an aquifer, the electrical resistivity measurements are influenced by the fluid contained within the aquifer which adds a crucial aspect to the analysis.Within the study area, a total of 30 evenly distributed VES probes were conducted using the DC Terrameter SAS 4000 instrument of ABEM, Sweden.The instrument allows us to measure the apparent resistivity of subsurface strata.The locations of VES probes were recorded using a GPS.The survey involves passing a low frequency electrical current through two current electrodes (A and B) and measuring the potential difference between two potential electrodes (M and N) firmly fixed in the ground (Fig 5a).The VES technique is described as keeping the central point of the electrode array fixed while the distance between electrodes is increased systematically to gather information from deeper subsurface sections.Throughout the survey, all four electrodes are kept collinear while the potential electrodes are always fixed within the current electrodes.
In a uniform subsurface, the depth of electric current penetration is directly related to the electrode distance, and altering the electrode distances provides a more insightful view of subsurface stratification [48,49].Different electrode configurations can be employed depending on the investigation type, field conditions, and the sensitivity of the instrument.Although both Wenner and Schlumberger electrode configurations are commonly used, however, the Schlumberger configuration is more suitable for groundwater investigations [50][51][52].Thus, the Schlumberger electrode configuration (Fig 5a) was used for the study area with a current electrode separation (AB/2) ranging from a minimum of 1.5 m to a maximum of 300 m.The  weakening signals of the resistivity meter with each increase in AB/2 spacing.Consequently, the MN/2 spacing was increased, allowing for two readings of the same AB/2, one with a short MN/2 spacing and the other with a long MN/2 spacing.The potential electrode separations (MN/2) were established as 0.5, 5, 15, 30, and then 50 m.Expanding MN/2 enables readings to be taken using the same current electrode spacing (AB/2) for both the previous and expanded MN/2 values.
The collected VES data were then analyzed using IPI2Win computer software [8,5,[53][54][55].This software incorporates the collected field data to generate the resistivity model (Fig 5b).The model is created by comparing the field data with synthetic data derived from the processing, aiming to minimize the root mean square (RMS) error [3,56,57].The iteration process continues until the fitting errors between the field and synthetic data reach a minimum and remain consistent.Fig 6 presents a summary of VES data interpretation conducted in the study area along with the corresponding RMS errors.The RMS errors of all VES survey data range between 0.32% and 1.31% which is considered satisfactory (Fig 6).In order to create two-dimensional (2D) models depicting the resistivity values at various depths, VES points are acquired in nearby boreholes or functioning tube wells.These models consist of distinct horizontal layers with well-defined resistivity bands.
We have conducted a thorough examination of the RMS errors for VES models (Fig 7).The RMS errors for 29 out of 30 VES models range from 0.32% to 0.96%, showing exceptional accuracy, with the majority below 1% (Fig 7).Notably, VES-30, while displaying accuracy, exhibits a slightly higher RMS error of 1.31% (Fig 7), albeit still within reasonable bounds.This underscores the robustness and reliability of the derived resistivity models, affirming their suitability for subsurface characterization.
Each model of a processed VES probe represents a hypothetical VES response of horizontally stratified earth layers, typically comprising a limited number of sediment layers.It is worth noting that the base layer in each model extends to an unspecified depth allowing for flexibility in interpreting the subsurface structure.This information was utilized to create 2D and 3D maps that illustrate the variations in electrical resistivity at different depths across the study area.A state-of-the-art computer software Arc GIS version 10.5 was employed to generate resistivity contours and 2D maps of the study area.
Dar-Zarrouk (D-Z) parameters are indispensable for understanding the underground hydrogeological characteristics and identifying zones with aquifer potential.These crucial parameters, comprising longitudinal conductance (Sc), longitudinal resistivity (ρL), and transverse resistance (Tr), are computed using equations.
Sc represents longitudinal conductance measured in mho, ρL stands for resistivity in Ω-m, and Tr denotes transverse resistance in Ω-m 2 .To calculate these values, ρ represents the average resistivity of saturated layers in Ω-m, while h signifies the average thickness of subsurface saturated layers in meters.These numerical values are derived from VES models.The distribution of Dar-Zarrouk parameters across a geographical area is visualized using ArcGIS 10.5, employing the interpolation method.

Physiochemical analysis
To provide a comprehensive overview of the most recent groundwater conditions of the study area, a direct investigation employing physiochemical analysis of groundwater samples was also conducted in addition to the ERS method of indirect investigation.In order to collect groundwater samples for field testing and further detailed physiochemical assessments in the lab, a field survey in the research area was carried out in October 2022.A total of 20 groundwater samples were initially collected and dispersed across the widespread region at depths between 30 and 180 feet (Fig 1).The TDS and pH levels were instantly examined in the field using standard portable equipment.The groundwater samples were also adequately preserved using Polyethylene terephthalate bottles for laboratory physiochemical analysis.To measure the concentrations of essential anions and cations, ion chromatography, ultraviolet-visible spectrophotometer, and several standard analytical methods were used.The graphite furnace atomic absorption method was used to calculate the concentration of arsenic present.Arsenic was measured in ppb or μg/l while the TDS and other physiochemical parameters were recorded in mg/l.The samples were chosen for comprehensive analysis in inhabited regions, specifically to assess the latest quality of drinking water for the residents.The physiochemical analysis was done in accordance with the WHO standards for different parameters to check the fitness of groundwater for drinking usage.

Interpretation and modeling of VES data
In order to gain a comprehensive understanding of hydrogeological layers in the study area, a thorough analysis was conducted on the collected VES data.This analysis revealed the presence of distinct four to six subsurface layers, each of them characterized by unique resistivity values with depth/thickness profiles.These resistivity measurements played a vital role in identifying key lithological properties such as particle size, porosity, and fluid content of subsurface strata.To establish a reliable correlation between lithology and resistivity, the processed data of 30 VES probes were meticulously examined alongside lithological logs obtained from nearby boreholes.Table 1 summarizes this valuable information while Fig 8 provides a visual representation of careful examination enabling the calibration of resistivity values and validation of the lithological interpretations.The 2D/3D models were then developed by analyzing all VES curves and using the information from Table 1.Based on previous literature, borehole information, reconnaissance field visits, surface geology, and hydrological setup of the research area, the resistivity values are classified as very low to very high (Table 1).
In a nutshell, the subsurface alluvial deposits in the study area are comprised of a diverse range of sediment types including intermixed clay, silt, sand, gravel, and some kankar which were encountered at various depths within aquifer-bearing boreholes (Fig 8).It is also worth noting that the resistivity inversion results confirm the existence of a localized region with clay-rich sediments resulting from salt intrusion, indicating the presence of aquifer zones with deteriorating water quality below the water table.However, differentiating resistivity readings between shallow and deep alluvium posed challenges during the interpretation of the VES results due to the relatively low resistivity contrast.Nevertheless, the borehole data (Fig 3) provided supporting evidence for the predominance of sand-gravel sediments as the primary aquifer materials for fresh/good water resources.
The analysis of the digital elevation map (Fig 4) revealed a gradual decline in elevation from northwest to southeast which corresponded to the flow direction of streams.This elevation pattern exhibited a correlation with variations in groundwater levels across the study area.For instance, borehole TH-8 in the eastern region indicated a water table depth of approximately � 5m while TH-7 in the northwest indicated a depth of around � 3.5m.In areas where the penetration depth of the VES was limited, the resistivity data indicated a higher prevalence of clay-rich sediments, impeding the flow of electric current at greater depths.Dry sand-gravel deposits showed elevated resistivity readings at the surface and above the water table while saturated sand-gravel sediments showed high resistivity readings below the water table.Predominantly sand-gravel sediments saturated with water comprise the majority of the primary aquifer of fresh/good quality water in the study area with resistivity values 40 ->230 O.m.The effective management and exploitation of groundwater assets can  be improved significantly by these findings, providing insightful information on underlying lithology and groundwater conditions.A comprehensive interpretation was carried out to ascertain the true resistivity measurements and thicknesses of distinct subsurface layers in order to gain an in-depth understanding of resistivity distribution in the study area.Following a systematic interpolation and visual representation of these results on a base map, the 2D resistivity maps at various depths (2 m, 10 m, 25 m, 50 m, 75 m, 100 m, 125 m, 150 m, 175 m, 200 m, 225 m, 250 m, 275 m, and 300 m were generated (Fig 9).These maps played a vital role in assessing the subsurface groundwater conditions.The resistivity values recorded in the study area exhibited a wide range, varying from 11.5 to 634.9 Ω-m simultaneously.The depth of the water table showed fluctuations between 3.5 and 6 m across different locations.
To categorize these resistivity values based on their correlation with water presence, the upper strata above the water table specifically at a depth of 2 m were thoroughly examined (Fig 9a).At this depth, very low resistivity values (less than 20 Ω-m) are the indication of fine silts/ clays.As the conditions are above the water table, so, the low resistivity values (20 -<40 Ω-m) are associated with the presence of dry sediments primarily composed of sand with some silt and clay.Medium resistivity values (40 -<100 Ω-m) are associated with dry sediments predominantly comprised of medium to coarse sand.Sediments exhibiting high resistivity values ranging from 100 Ω-m to 230 Ω-m mainly consist of dry sand with a mixing of thin layers of silt/clay.Furthermore, resistivity values surpassing 230 Ω-m, indicate very high resistivity zones characterized by predominantly coarser sand and gravel.
At shallow and greater AB/2 values, ranging from 10-300 m (Fig 9b -9n), the focus is shifted to investigate the hydrogeological conditions below the water table.At these depths, the resistivity values below 20 Ω-m are interpreted as the presence of sediments primarily composed of silty clay or clayey silt indicating saline/brackish groundwater saturation.Sediments with low resistivity values (<40-20 Ω-m) were interpreted as sandy layers with thin layers of silts and clays implying the saturation of marginally suitable quality of groundwater.Within the resistivity range of <100 to 40 Ω-m (medium resistivity) at these depths, the sediments are predominantly characterized by medium to coarse sandy layers with minor gravel indicating good/fresh quality groundwater saturation.The high resistivity range (<230-100 Ω-m) is mainly associated with sandy layers containing gravel indicating good quality groundwater saturation.Minor or rare intermixing of fine materials such as silts and clays may also be associated with the strata of high resistivity.Resistivity values exceeding 230 Ω-m are interpreted as the existence of sand, gravel, and stiff kankar inclusions, signifying again the conditions of good quality groundwater saturation.However, the stiff kankar inclusions may compromise the water-well yield to some extent at some subsurface levels depending upon the stiffness of the layers.To gain a comprehensive understanding of hydrostratigraphy in the study area, we focused on a specific profile labeled A-A' on the base map ( Fig 10).predominantly fresh/good from shallow depths of about 8m down to the maximum investigation depth (300m), with only minor localized groundwater occurrences of moderately fair and brackish quality.

By interpolating the resistivity data along profile
The groundwater recharge in the study area primarily relies on the River Indus which flows parallel to the western region.Canals and to some extent the precipitation also contribute to the recharge process.The western side of the study area may get the benefits from a robust groundwater recharge due to the proximity of the River Indus while the recharge gradually decreases as we move away from the river towards the east.Nevertheless, there are some significant pockets with groundwater of moderately acceptable or marginally suitable quality extending up to the depth of 300 m (Figs 9 & 12).To maintain a sustainable water balance and preserve the quality of fresh groundwater zones, it is better to utilize groundwater resources up to the shallow depth of about 50 m within the identified areas of moderately acceptable groundwater quality.It is essential to avoid water extraction from saline/brackish zones below the depth of 50 m, because such activity may produce a risk of disrupting water balance due to limited recharge at deeper levels.Notably, a significant portion of the study area beyond a depth of 150 m (Fig 9a -9n), exhibits fresh/good quality groundwater making it suitable for extraction purposes.
Dar-Zarrouk parameters (longitudinal conductance, longitudinal resistivity, and transverse resistance) maps are produced through the interpolation of numerical values extracted from

Physiochemical modeling
Evaluating the quality of groundwater is of utmost importance to assess its suitability for various purposes, such as drinking, domestic use, agriculture, and industry [58].In order to thoroughly investigate this matter, we conducted a comprehensive analysis by collecting water samples from 20 different locations within the research area (Fig 1).These samples were obtained at depths ranging from 30 to 180 ft with an average depth of 69 ft.To determine the physicochemical characteristics of groundwater, the samples were carefully analyzed in the laboratory.The results of this analysis were then compared against the water quality guidelines established by the WHO [59][60][61] (Table 2).7.9 with an average value of 7.4 (Table 2, Fig 14b).The moderately alkaline nature of groundwater across all climate conditions can be attributed to factors such as CO 2 loss and accumulation of mineral salts.The analysis of total dissolved solids (TDS) revealed a range of 465 to 3817 mg/l with an average concentration of 989 mg/l (Table 2).Notably, approximately 30% of samples exceeded the maximum acceptable limit of 1000 mg/l indicating unsuitability for drinking purposes (Fig 14c).In terms of arsenic content, the average concentration in the groundwater of the study area was found 18 μg/l which exceeds the permissible limit (Table 2).It is worth mentioning that seven out of 20 samples exhibited significantly higher levels of arsenic than the allowable limit (Fig 14d).The average fluoride content in the groundwater was 1.1 mg/l with a range of 0.2 to 2 mg/l (Table 2).Approximately 20% of the water samples exceeded the WHO guidelines for fluoride (Fig 14e).Similarly, the average iron content in the groundwater samples was 0.3 mg/l ranging from 0.1 to 0.7 mg/l (Table 2).With regard to the water quality standards for Fe 2+ levels, seven out of the 20 sites exceeded the guidelines set by WHO (Fig 14f).The concentration of nitrate (NO 3-) in the groundwater samples ranged from 10 to 90 mg/l with an average of 34.5 mg/l (Table 2).Normally, natural groundwater contains nitrate levels below 10 mg/l but elevated proportions were observed in certain zones (Fig 14g ) indicating the influence of nitrogen-rich fertilizers and human wastes as well.The concentration of nitrite (NO 2

−
) in the groundwater samples varied from 0.1 to 1.0 mg/l (Table 2) and about 15% of the water samples exceeded the WHO guidelines for nitrite concentration in drinking water (Fig 14h).Human civilization and the usage of agricultural fertilizers in the region are responsible for the high nitrite concentrations in some samples.
A common and straightforward approach for quantifying the linear relationship between two variables is the utilization of the Pearson correlation coefficient.This coefficient, which ranges from -1 to +1, encapsulates possibilities of positive, negligible, or negative correlations.We employ it to comprehend both the strength and direction of associations within our data.Within our investigation, we have computed Pearson correlation coefficients to assess the interplay among multiple water quality parameters encompassing depth, TDS, pH arsenic, iron, fluoride, nitrate, and nitrite, as illustrated in Fig 15 .These correlation coefficients unveil potential links between diverse facets of water quality Fig 15.This analysis serves as a crucial step in identifying potential sources of contamination, natural fluctuations, and the underlying mechanisms governing these variables within aquatic ecosystems.The robust correlations observed among various parameters are attributed to shared origins, hydrogeological and geomorphic characteristics, human activities, and intricate chemical interactions.Natural geomorphic conditions, mineral dissolution, and weathering processes contribute to the release of ions into water bodies, thereby affecting their concentrations.Simultaneously, anthropogenic inputs such as domestic wastewater and fertilizers introduce pollutants into water sources in the study area, exerting an influence on the observed correlation patterns.Additionally, the  complex chemical interactions between ions contribute to equilibrium shifts and concentration variations, further contributing to the observed correlation patterns.
The Schoeller diagram analysis continuously exhibited a clear pattern in the water specimen curve, illustrating a high degree of consistency across the groundwater sources in the study area (Fig 16).Anthropogenic influences and regional hydrological factors are responsible for the consistency in pH, TDS, arsenic, fluoride, nitrite, nitrate, and iron concentrations found in groundwater supplies within the study area.The use of groundwater resources within the area must be sustainable, therefore these results have practical implications for detecting possible sources of contamination, putting into practice the efficient water management plans, and ensuring these measures are implemented.
In order to better comprehend the relationships between various physical and chemical variables, we employed a scatter plot matrix to graphically analyze the correlation, strength, and orientation among different parameters (Fig 17).These findings suggest that the physiochemical variables under investigation are not interrelated or have a significant impact on each other.The values of each parameter appear to be impacted by multiple sources or mechanisms, leading to independent deviations among sample locations.This may imply that several causes or sources are contributing to the contamination levels.

Conclusions
The findings of this research work yield valuable insights into the subsurface hydrogeological layers and groundwater conditions in the study area.These insights were derived through the interpretation of VES data and physiochemical analysis, shedding light on the distinct subsurface layers characterized by unique resistivity values and depth/thickness profiles.Four to six geo-electric subsurface layers were identified through the processing of VES data.The The spatial distribution of resistivity was analyzed by generating 2D resistivity maps at various depths allowing the identification of fresh/good, marginally suitable, and saline/brackish groundwater zones.Further insights into hydrostratigraphy were gained by examining resistivity values and lithological cross-sections along a selected profile.Groundwater recharge primarily occurs through the River Indus with more significant recharge observed near the river.
Physiochemical analysis of groundwater samples indicated slightly moderate alkalinity in some parts of the study area.Moreover, the analysis revealed varying levels of TDS, occasional samples surpassing permissible limits for arsenic contamination, elevated concentrations of iron and fluoride in some specific locations, and some parts with higher content of nitrate and nitrite levels influenced by human activities and fertilizer use.The analysis identified the areas where groundwater is unsuitable for drinking due to elevated TDS, arsenic, and fluoride levels exceeding the recommended guidelines.Based on these findings, it is advisable to utilize groundwater resources up to a depth of about 50 m in the zones with good quality groundwater while refraining from extraction in brackish zones.Extensive areas containing good quality groundwater suitable for extraction are also present at greater depths.Ultimately, this research contributes vital knowledge on subsurface lithology and groundwater conditions facilitating the effective management and utilization of groundwater resources in the study area.
Fig 1 represents the location of VES probes, boreholes drilled in the study area, and collected groundwater samples.

Fig 6 .
Fig 6. (a-f) Modeled resistivity curves.Apparent resistivity data are marked by small black circles.The red curve is the best-fitted line with apparent resistivity data.The solid blue block line is the modeled resistivity (synthetic resistivity).The horizontal axis is the current electrode spacing (AB/2) in meters and the vertical axis is resistivity in ohm meter.https://doi.org/10.1371/journal.pone.0302442.g006

Fig 8 .
Fig 8. Calibration between lithology and resistivity data demonstrated through the visualization of borehole litho-logs and modeled resistivity curves.https://doi.org/10.1371/journal.pone.0302442.g008 A-A' (Fig 10), we have developed a resistivity cross-section referred to as A-A' (Fig 11a).Additionally, Fig 11b provides an interpreted lithological cross-section, offering a detailed two-dimensional representation of lithological variations along this profile.The generation of a 2D cross-section (Fig 11a) along the designated profile involved utilizing processed electrical resistivity data in conjunction with the closest borehole data to profile A-A'.The resistivity values along profile A-A' exhibit significant fluctuations encompassing both lower and higher values.Above the water table, the resistivity values displayed spatial variations from the western to the eastern side of the profile (Fig 11a).The upper section primarily consisted of dry sediments such as sand, silt, and clays.Below the water table, bluish colors (light and dark blue regions) in Fig 11b correspond to yellow and orange colors in Fig 11a, indicating the water saturated

Fig 12 .
Fig 12. Enhanced visualization of true resistivity distribution (An aerial view from SE &SW): Stacked multi-layered maps organized according to resistivity ranges.https://doi.org/10.1371/journal.pone.0302442.g012 Fig 14 represents the 2D visualization of the results of different physiochemical parameters which further lead to the differentiation of the quality zones of groundwater in the study area.The pH level of groundwater samples in the study area exhibited variations between 7.1 and